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Molecular dynamics simulations are used to investigate the effects of deformation on the segmen- 
tal dynamics in an aging polymer glass. Individual particle trajectories are decomposed into a series 
of discontinuous hops, from which we obtain the full distribution of relaxation times and displace- 
ments under three deformation protocols: step stress (creep), step strain, and constant strain rate 
deformation. As in experiments, the dynamics can be accelerated by several orders of magnitude 
during deformation, and the history dependence is entirely erased during yield (mechanical reju- 
venation). Aging can be explained as a result of the long tails in the relaxation time distribution 
of the glass, and similarly, mechanical rejuvenation is understood through the observed narrowing 
of this distribution during yield. Although the relaxation time distributions under deformation are 
highly protocol specific, in each case they may be described by a universal acceleration factor that 
depends only on the strain. 
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I. INTRODUCTION 

Although it has long been recognized that mechani- 
cal deformation can modify the structural relaxation dy- 
namics of polymer glasses, obtaining a complete picture 
of the deformation-induced mobility remains an impor- 
tant challenge. It is most often assumed that mechanical 
stress lowers the activation barriers to molecular mobil- 
ity, thereby allowing the solid to yield and flow. This 
simple idea, first proposed by Eyring in 1936 [l[ , remains 
the central tenet of several models of plasticity in amor- 
phous solids (HQ. Other control variables have also been 
proposed to explain the increased molecular mobility un- 
der deformation. In the Soft Glassy Rheology (SGR) 
model, activated barrier crossings are accelerated by the 
local strain on mesoscopic regions of the glass 0, Q . In 
the Shear Transformation Zone (STZ) theory, deforma- 
tion is caused by a population of bistable states (shear 
transformation zones) that can flip from one orientation 
to another in response to stress @, 0). The stress in- 
duced transition between the bistable states is modeled 
using a variant of the Eyring model, and the population 
of shear transformation zones is controlled by an effec- 
tive temperature which explicitly depends on the strain 
rate e [a, M . Still other models have proposed the free 
volume jlO|. and the configurational entropy (TT| as the 
important variables to describe the deformation-induced 
mobility. There is clearly no consensus on the correct 
way to think about deformation in glasses. 

Experimentally, the structural relaxation times of de- 
formed polymer glasses have traditionally been studied 
through rheological measurements, where the segmen- 
tal dynamics are inferred from the mechanical response 
[Tol [T2I - FT5I . In recent years, experiments and simula- 
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tions have allowed researchers to directly probe the mi- 
croscopic segmental dynamics during deformation. Solid- 
state NMR experiments examined the dynamics of glassy 
Nylon-6 near T g under uniaxial extension and found 
significant deformation- induced mobility near yield [l6j . 
Lee et al. performed a series of experiments investigating 
the segmental relaxation dynamics in solid poly(methyl 
methacrylate) through an optical photobleaching tech- 
nique that probes the reorientation of small dye molecules 
embedded in the polymer matrix during a step stress 
(creep) experiment [171 . 1 1 81 ] . Relaxation times were found 
to be significantly accelerated compared to the unde- 
formed sample, and in fact the segmental mobility was 
close to linearly correlated with the strain rate. During 
plastic flow, the distribution of relaxation times was also 
narrowed, which was interpreted as evidence that the dy- 
namics has become more homogeneous. 

Molecular dynamics simulations have confirmed these 
results through direct measurement of the segmental dy- 
namics. Lyulin et al. performed simulations on united 
atom models for polystyrene and pol ycarbonate during 
constant strain rate deformation jl9j ]. The segmental 
mobility was monitored through the monomer mean- 
squared displacement and increased dramatically beyond 
the yield point. Similarly, simulations of polyethylene 
under compressive strain showed that the dihedral tran- 
sition rate in between trans/gauche states was increased 
due to the deformation even in the sub-yield regime [20| . 
Several simulation studies have also confirmed the strong 
correlation between the relaxation time and strain rate 
that has been seen in experiment. Riggleman et al. simu- 
lated a bead-spring polymer glass during non-linear creep 
[2ll |22| . The covalent bond autocorrelation and interme- 
diate scattering functions both decayed more rapidly in 
the deformed sample, and it was found that the segmen- 
tal mobility was strongly correlated to the strain rate 
under both tensile and compressive loading. 

For polymer glasses in particular, the effects of physi- 
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cal aging further complicate the analysis of the dynamics 
under deformation, as the mechanical response depends 
sensitively on the history of the glass. Non-equilibrium 
relaxations cause relaxation times to grow, and hence the 
glass becomes stiffer with the time since it was formed. 
Since the age of a sample is often characterized by its 
relaxation times, the deformation-induced acceleration 
of the dynamics has been associated with "mechanical 
rejuvenation" of aging glasses [Tof. Indeed, in the post- 
yield flow regime, all history dependence is erased. Hasan 
and Boyce investigated this phenomenon using differen- 
tial scanning calorimetry on polymer glasses during con- 
stant strain rate deformation |23|. While, at sub-yield 
strains, the aged glasses exhibited different mechanical 
properties and enthalpies, these differences were found 
to disappear once the glasses began to flow. Simulations 
have provided further insight into the nature of rejuvena- 
tion. Utz et al. found that the history dependence of the 
inherent structure energies and the pair correlation func- 
tions of the binary metallic Lennard- Jones glass (BMLJ) 
were also erased during yielding [24j . Post-yield, these 
parameters took values more typical of rapidly quenched, 
or unaged samples. Further analysis, however, seems to 
indicate that the rejuvenated glass is not truly the same 
as a younger glass |25|, [26| . 

Below yield, the claim of mechanical rejuvenation is 
even more controversial 1271 ] . In this regime, deforma- 



tion increases molecular mobility at the same time as 
aging decreases it [HI, and experiments show that both 
the volume [2!| and the relaxation times [28| of the glass 
quickly return to their previous states once the load is 
removed (30j . In this report, we reserve the term "re- 
juvenation" for the erasure of aging due to deformation, 
without comment on whether the yielded state truly re- 
sembles a younger glass or is something else altogether. 

Mobility and aging in previous studies has been mea- 
sured using spatially averaged correlation and response 
functions; however, it is becoming increasingly clear that 
glasses are both dynamically and mechanically heteroge- 
neous. A number of simulations have studied the spatial 
heterogeneity of structural relaxations under deformation 
in the limit of low temperatures and small shear rates 
[3ll |32| . Results show that the elementary relaxations 
responsible for plastic deformation are strongly localized, 
involving a region of several diameters across. They are 
also intermittent: rapid relaxation events are separated 
by periods of perfectly affine deformation. Furthermore, 
the stress and elastic modulus are also locally heteroge- 
neous [33L [111 and structural relaxations occur predomi- 
nantly in weak regions of the glass [33[ . 

We use molecular dynamics simulations of an ag- 
ing polymer glass to investigate the complex interac- 
tion between aging and deformation at the microscopic 
level, and to identify the deformation variables impor- 
tant to the description of the deformation-induced mo- 
bility. Three deformation protocols are employed: a step 
stress (creep), a constant strain rate deformation, and 
a step strain (stress relaxation). The microscopic relax- 



ation events are identified as discontinuous "hops" in the 
individual monomer trajectories which occur in between 
periods of affine deformation, and the full distribution 
of hop times and displacements are evaluated. Previ- 
ously, we reported that the changes to these distributions 
in three different deformation protocols and for different 
ages exhibit a universal dependence on global strain [35j . 
In the present paper, we provide a more detailed account 
of these findings, present additional data, and expand 
on the analysis. In Sections [TT] and IIIIl we outline the 
details of our simulations and our methodology for de- 
tecting hops. Section IIVI presents the results of these 
simulations, and in Section [V] we discuss our finding and 
conclusions. We also present in the Appendix further de- 
tails of the hop analysis technique used to determine the 
relaxation time distributions. 



II. SIMULATIONS 

We perform molecular dynamics (MD) simulations of 
a well-known bead-spring polymer model for polymer 
glasses [H, H3|- The beads interact via a non-specific 
van der Waals interaction given by a 6-12 Lennard- Jones 
potential, and the covalent bonds are modeled with a stiff 
spring that prevents chain crossing (38|. The Lennard- 
Jones interaction is truncated at 1.5 times the bead di- 
ameter and adjusted vertically for continuity. Unless oth- 
erwise noted, all polymers have a length of 10 beads, and 
each sample contains two thousand polymers. To check 
that results are insensitive to chain length, we also ex- 
amined selected cases with chains of length 100 (slightly 
longer than the entanglement length [331). The simula- 
tion box is initially cubic with periodic boundary condi- 
tions. All results will be given in units of the diameter a 
of the bead, the mass m, and the Lennard- Jones energy 
scale, 1*0. The reference time scale is tlj = yjma 2 /u$. 

To create a glass, the ensemble of polymers is first re- 
laxed at constant volume and a melt temperature of 1.2 
and then cooled to form a solid at T — 0.25 (T g k, 0.35 
40]). The density of the melt is chosen such that af- 
ter cooling the pressure on the box is zero. The glass 
is then aged for various wait times t w at constant tem- 
perature and zero pressure using a Nose-Hoover thermo- 
stat /barostat. For the creep protocol, a constant uniax- 
ial tensile stress in the x-direction is ramped up over a 
time of 75 and then held constant using the Nose- Hoover 
barostat. Four different stresses are applied a = 0.2, 0.3, 
0.4, and 0.5 in order to investigate a range of behaviour 
from the linear regime to yielding and flow. In the sec- 
ond protocol, a constant shear rate is applied by rescaling 
the particle positions and the simulation box in the x di- 
rection. Three strain rates are used: e = 2.2 x 10~ 7 , 
8.9 x 10~ 7 , and 8.9 x 10~ 6 . For the final protocol, a step 
strain is applied over a time of 75 and then held constant. 
Strains of 0.01 to 0.04 are evaluated. In all three cases, 
the stress is maintained at zero in the other two direc- 
tions, and their dimensions respond dynamically to the 
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FIG. 1. Non-affine displacements Ax, Aj/, and Az 
(black/blue/red) for a monomer undergoing creep at a — 0.5. 
The vertical dotted lines indicate relaxations identified by our 
hop detection algorithm. The total strain during this time is 
approximately 0.8. 



the displacement is often much larger than the cage size. 
The ambiguity in determining a hop from a single trajec- 
tory is related to the fact that the single particle hop is 
only a small part of a collective relaxation event. In this 
work, we therefore identify hops through a threshold in 
the standard deviation of the displacement, rather than 
the mean; that is, we define relaxation events by the ac- 
tivity of the particle, rather than the final state. Further 
details of the procedure are discussed in the Appendix. 
Once the trajectories are reduced to a series of hops, we 
construct the probability distribution functions for the 
first hop time p{t\) (the time from the beginning of the 
measurement to the first hop), the persistence time p(r) 
(the time in between all subsequent hops or the residence 
time for a caged configuration), and the hop displace- 
ments p(dx). 



IV. RESULTS 
A. Aging at zero pressure 



imposed deformation through the Nose-Hoover barostat. 



III. HOPPING DYNAMICS 

In addition to monitoring the mechanical response in 
the deformation protocols described above, we are inter- 
ested in understanding more directly how the deforma- 
tion accelerates the particle relaxation dynamics. Struc- 
tural relaxations are identified on the scale of a single 
particle or monomer segment by monitoring a subset of 
five thousand particle trajectories chosen randomly from 
the simulation volume. Figure [1] shows an example of 
the non-affine particle displacement during a creep ex- 
periment (the portion of the displacement due to the 
macroscopic deformation of the simulation volume has 
been removed). The trajectories take the form of long 
periods of vibration where their motion is predominantly 
affine, punctuated by discrete hopping events. 

Single particle hopping dynamics have been studied in 
simulations of metallic liquids and glasses [ill . l42j , in ki- 
netically constrained models of supercooled liquids (43l |. 
and in experiments of colloidal glasses 44| and granu- 
lar materials [45l j . and have been very useful in under- 
standing the nature of heterogeneous dynamics in soft 
glassy materials. To detect the individual hops, a run- 
ning average and standard deviation is computed for the 
non-affine displacement of each particle over a time of 
400. The usual method of defining a hop is through a 
threshold in the displacement of the mean particle po- 
sition, where the threshold is determined in relation to 
the vibrational amplitude [iH I43H45} . However, we have 
observed that frequently particles experience large, inter- 
mittent fluctuations in position with little net displace- 
ment. These events can not be disregarded as hops, as 



At finite temperature, the undeformed glass continues 
to evolve in a non-ergodic manner due to activated ther- 
mal processes. The non-equilibrium relaxation dynamics 
causes aging in all of the properties of the glass with the 
wait time t w since vitrification. The volume decreases 
logarithmically with the wait time, as does the potential 
energy. The most striking results of aging are, however, 
seen in the dynamical properties of the glass. Glasses 
obey "time- waiting time superposition" : two-time corre- 
lation and response functions can be rescaled by the wait 
time C(t, t w ) = C vib (t) + C con f(t/t^), where fi is the ag- 
ing exponent and t is the measurement time. Although 
glasses are characterized by a wide distribution of relax- 
ation times, it appears as though they are all rescaled 
identically during aging. 

Previously, we reported on the hop statistics in a poly- 
mer and binary metallic Lennard- Jones glass undergoing 
simple aging after a rapid quench [46] . Hop statistics for 
the polymer glass can be seen here again as solid lines in 
Figs, gfb-d), and also in Figs. ISfb-d) and EJb-d). The 
first hop time distributions (Fig. [3](b) ) take the form of 
two power laws, and the intersection between the two in- 
creases with t$j in the same way as the correlation and 
response functions. This time seems to be correlated with 
the cage escape time r Q where structural relaxations be- 
come important in the glass. Alternatively, Figure OHc) 
shows that the distribution of persistence times decays 
as a single, age- independent power law, p(r) ~ t -12 . 
In ref. [46j, we showed that aging can be understood as 
the direct result of the wide distribution of persistence 
times which, in the thermodynamic limit, has an infinite 
mean. Evolving the system of particles using a continu- 
ous time random walk parameterized with the initial con- 
dition p(ti,t w = 0), which reflects the quench, and the 
persistence time distribution p(r) is sufficient to predict 
the observed changes to the first hop time distributions 
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FIG. 2. Mean-squared displacement as a function of the num- 
ber of hops n that a particle has experienced for the polymer 
model (o) and the BMLJ model (□) (simulation details of 
the BMLJ results can be found in [4(j]). Lines have slopes of 
1/2 (blue) and 1 (red). Inset shows the corresponding mean- 
squared displacement for both models with time. 



p(h,t w ) during aging. 

A power law form for p(r) was also found for the 
BMLJ supercooled liquid by using a rigorous energy land- 
scape approach to detect transitions [47j- Doliwa and 
Heuer performed frequent energy minimizations in or- 
der to detect collective relaxations in between adjacent 
metabasins in the potential energy landscape. The per- 
sistence time distributions were found to take the form 
of two power laws: at short times the transitions are con- 
tinuous and liquid-like, and at long times the transitions 
are dominated by uncorrelated hops between metabasins. 
The decay at long time becomes less steep with tempera- 
ture, approaching t~ 2 as the glass transition temperature 
is approached from above; that is, (t) is finite in the liq- 
uid state. Vollmayer-Lee measured the persistence time 
distributions in the binary Lennard- Jones glass using a 
single particle hop detection algorithm which differenti- 
ates between hops that lead to an entirely new position 
(irreversible hops) and hops back and forth between two 
average positions (reversible hops) [4ll ]. p{r) was again 
found to decay in a weak power law (t~ 0,84 ), independent 
of history and temperature in the glassy state. 

Figure [3Jd) shows that the hop displacement distri- 
bution p(dx) also does not depend on the wait time. 
p{dx) is sharply peaked around zero and decays rapidly 
at large dx. The rapid decay of the tails results from 
the connectivity of the polymer glass; the tails of p(dx) 
are purely exponential in the binary metallic Lennard- 
Jones (BMLJ) glass (4(|. Polymer specific dynamics is 
also apparent in the diffusion behaviour of the polymer 
glass. The mean-squared displacement (Ax(t) 2 ) of both 
the polymer model and the binary metallic Lennard- 



Jones glass show the two-step relaxation behaviour char- 
acteristic of glassy materials (inset Fig. [5]). However, in 
the cage-escape regime, (Ax(t) 2 ) increases much more 
rapidly in the metallic glass than in the polymer glass. 
The reason for this is evident if the mean-squared dis- 
placement is plotted versus the number of hops n a par- 
ticle has experienced rather than the time elapsed t, as 
shown in Figure O The mean-squared displacement of 
the polymer beads increases with the number of hops as 
(Ax(n) 2 ) ~ ^fn (Rouse dynamics) rather than ~ n for 
simple molecular glass formers due to the constraints of 
the covalent bonds. Alternatively, the aging behaviour of 
the relaxation time distributions were found to be qual- 
itatively similar for the polymer model and the BMLJ 
glass [4(3]. This is a good indication that the results dis- 
cussed here apply quite generally to glassy materials, and 
the polymer-specific aspects of the dynamics are seen pri- 
marily in the displacements rather than the relaxation 
times. 



B. Step stress 

We next turn our attention to the effects of deforma- 
tion on the relaxation dynamics. Figure [3] shows the me- 
chanical response (a), and the hop statistics (b-d) in the 
step stress experiment. When a stress a is suddenly ap- 
plied to the glass, the strain e(i, t w ) shows an initial elas- 
tic response followed by slow elongation, or creep, due 
to structural relaxations in the glass. In Figure EJa), 
the creep compliance J(t,t w ) — e(t,t w )/a is plotted as 
a function of wait time and stress. At small stress, the 
creep compliance takes the form of a stretched exponen- 
tial at short times, changing to a logarithmic increase at 
t ^> t w due to the effects of aging at long times [13, 51] ■ 
For increasing wait time, the glass becomes stiffer and 
the creep compliance curves at short times can be su- 
perimposed to form a master curve by rescaling the time 
variable with a power law in t w , 

J(t,t w )=J(t/t»). (1) 

The aging exponent /i is approximately 0.8 at a small 
stress of a = 0.3. At a = 0.5, Fig. [3^a) shows that the 
short time compliance depends on the wait time, but at 
longer times the glass yields and all history dependence 
is lost. The strain rate in this regime decreases over time 
due to the effects of strain-hardening. 

The hop statistics for the a = 0.4 and 0.5 creep ex- 
periments are shown in Fig. [HJb-d), and compared to the 
unstressed, aging glass. For a = 0.5, the first hop times 
are clearly accelerated by the stress. The tails of the dis- 
tribution decay with a much steeper power law than the 
undeformed glass it^ 2 ' 2 vs. £f ' ). This corresponds well 
to experiments, which showed that in the flow regime, 
not only are the segmental dynamics accelerated, but 
the width of the relaxation time distribution becomes 
narrower as well [Lsl ]. This distribution is narrowed for 
all stresses tested here, however the changes to p{t\) are 
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FIG. 3. Relaxation dynamics during the creep experiment, (a) Creep compliance for a — 0.4 (dashed), 0.5 (dotted). Squares 
show the a — 0.5 curve for chains of length 100. Distributions of (b) first hop times, (c) persistence times, and (d) displacements 
for the undeformed sample (solid lines), a = 0.4 (dashed), and a — 0.5 (dotted). In all plots t w = 750 (black), 7500 (red), and 
75000 (green). Circles in (c) show p(r) for t < T a ~ 10 4 (see text) and t w = 75000. Straight lines indicate power law with the 
given slopes. 



considerably smaller in the sub-yield regime. The per- 
sistence time distributions are also narrowed during the 
creep experiment. Unlike in the case of simple aging, p(r) 
during a creep experiment depends explicitly on both the 
wait time t w , and the measurement time t since the stress 
was applied. The distributions shown in Fig. E^c) are 
obtained from the persistence times of cages originating 
at times much greater than the cage escape time, where 
the strain is increasing logarithmically. The narrowing 
of p(r) is even more pronounced for cages originating at 
shorter times (circles). This result provides an interest- 
ing explanation for mechanical rejuvenation. Since aging 



can be understood as the effect of an infinite mean per- 
sistence time, then aging ceases when the tails of p{r) 
decay with an exponent steeper than —2, and a steady 
state is established which is independent of the age of the 
glass. For glasses below yield, the persistence time distri- 
bution remains in the aging regime, whereas at a = 0.5, 
a steady state flow regime has been established and the 
persistence time distribution decays with an exponent of 
2.2. 

Figure EJd) shows that the hop displacement distri- 
bution is wider for the a = 0.5 case. At this stress, 
p(dx) depends on the strain rate (discussed further be- 
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have a length of 10 beads, which is far shorter than the 
entanglement length for this model. It is important to 
check how the results depend on chain length in our poly- 
mer model, especially at large strains where strain hard- 
ening is pronounced in entangled polymer glasses. To 
this end, we have performed the same simulations with 
chains of length 100. Figure[3Ja) (squares) shows that the 
long-chain glass exhibits considerably more strain hard- 
ening during yield, with a lower strain and strain rate 
than the short-chain glass. We compare the results of 
the strain rate versus the hop rate for the two different 
chains lengths in Fig. @] Not only does the superposi- 
tion for different wait times and stresses still hold for 
the longer chains, but the form of the universal curve is 
identical to that found for the shorter chains. 



FIG. 4. Strain rate versus hop rate for a — 0.5 (o), 0.4 (□), 
0.3 (A), 0.2 (0) and t w = 750 (black), 7500 (red), and 75000 
(green). Solid symbols show results for chains of length 100. 
Inset shows the mean hop displacement versus the strain rate 
for the same data. 



low), and is highly anisotropic. The non-affine displace- 
ments in the direction of the applied stress are signifi- 
cantly smaller than displacements perpendicular to the 
stress. This does not appear to be true at any of the 
smaller stresses, and is likely an effect of strain hardening 
during flow. The fluctuations in the direction of strain 
may be suppressed in this regime because the chains are 
extended and under tension. 

Previous experimental [l?], [3 and simulation [2l|, [H[ 
studies have shown a strong correlation between the par- 
ticle mobility and the strain rate. The mobility was de- 
fined as the inverse of a relaxation time r obtained from 
the stretched exponential fits to an autocorrelation func- 
tion obtained periodically throughout the creep experi- 
ment, and was found to increase almost linearly with the 
strain rate. The analogous quantity to the mobility in our 
simulations is the average hop rate, which is determined 
by creating a histogram of the number of hops as a func- 
tion of time. A parametric plot of the strain rate versus 
the hop rate is shown in Fig. @] Data for all stresses and 
wait times fall on a universal curve, indicating that there 
is a direct relationship between the macroscopic strain 
rate and the microscopic hopping transitions. The strain 
rate does not increase linearly with the hop rate, however. 
Even at zero stress, the hop rate is finite due to thermally 
activated transitions. At higher strain rates, where the 
hopping dynamics are highly accelerated, the strain rate 
increases faster than linearly with the hop rate. This can 
be explained by the fact that the average hop displace- 
ment also increases with the strain rate, as shown in the 
inset. It is possible that this is due to unresolved hops at 
very short times, which is discussed further below. 

The polymer chains used to obtain the previous results 



C. Constant strain rate 

The stress-strain curve under constant strain rate de- 
formation is shown in Fig. [5][a). The effects of aging ap- 
pear as an increasing over-shoot stress, after which there 
is a period of shear softening and eventual flow. The size 
of this overshoot increases with the wait time, and with 
increasing strain rate [49| . As in experiments 23 1, the 
stress on our model glass in the flow regime is indepen- 
dent of the wait time; all history dependence has been 
erased. Because of the low strain rates studied here, the 
stress in the flow regime is also independent of the strain 
rate. 

The probability distributions describing the hop dy- 
namics for this protocol are shown in Fig. [5Jb-d). The 
constant strain rate deformation also greatly accelerates 
the hopping dynamics, but in a very different way than 
the step stress. At short times, the first hop time distri- 
butions are markedly unchanged, whereas, at long times, 
p{t\) no longer decays with power law tails but is in- 
stead sharply truncated at a timescale of approximately 
0.1 /e. The same behaviour is also evident in the persis- 
tence times. In contrast with the case of the step stress, 
during constant strain rate flow, p(r) does not depend on 
the wait time or the measurement time. This is interest- 
ing because in the flow regime strain hardening starts to 
become important, and yet the persistence time distribu- 
tion is the same at short times and at long times. Con- 
versely, the hop displacements are fairly isotropic at short 
times, and become wider and highly anisotropic during 
the flow regime. It would appear that strain hardening 
manifests in the spatial distribution of hops, rather than 
their relaxation times. 

Mechanical rejuvenation is seen after yield in the con- 
stant strain rate protocol as well, and can similarly be ex- 
plained by the persistence time distribution. The trunca- 
tion in p(r) at long times certainly means that the mean 
persistence time is finite and a steady state is established 
in the flow regime. 
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FIG. 5. Relaxation dynamics during the constant strain rate experiment, (a) Stress versus time for e = 8.9 x 10 7 (dashed) and 
e = 8.9 x 10 -6 (dotted). Distributions of (b) first hop times, (c) persistence times, and (d) displacements for the undeformed 
sample (solid lines) and e = 8.9 x 10~ 7 (dashed) and 8.9 x 10" 6 (dotted). In all plots t w = 750 (black), 7500 (red), and 75000 
(green). Straight lines indicate power law with the given slopes. 



D. Step strain 

The response to a step strain is defined by the modulus 
G(t,t w ) — a(t,t w )/e which is shown in Fig. (Ha). Upon 
applying the strain, there is an immediate increase in 
the stress of a glass, which is larger for glasses aged for 
longer wait times. As the strain is held constant, the 
stress gradually decays due to structural relaxations over 
a timescale which also depends on the wait time. Like 
the case of the creep compliance, the step strain modulus 
is shifted forward in time with increasing t w . For larger 
strains outside of the linear regime, G(t,t w ) decreases 
with increasing step strain e; the glass becomes softer 
because the dynamics are accelerated by the deformation. 

The first hop time distribution (Fig. [6jb)) shows that 
the number of relaxations occurring at very short times 
is greatly increased by the step strain; however, at longer 
times, p(ti) decays with the same power law as the un- 



deformed glass. At the largest strains investigated, p(t\) 
for all wait times appear to be almost identical, which 
would suggest a loss of history dependence (rejuvena- 
tion). The stress relaxation response indicates, however, 
that this is not the case. Even at very large strains, the 
modulus retains the features of aging. It is likely that 
the variation in p(t\,t w ) is to be found at short times 
where individual hops cannot be resolved by our meth- 
ods. Unlike the case of both the creep and the constant 
strain rate experiment, the persistence times (Fig. [BJc)) 
and the displacements (Fig. [^d)) are completely unaf- 
fected by the deformation: after a first hop, the particle 
loses all memory of the strain step. This result implies 
that the persistence times are not sensitive to the global 
stress. There is no time or wait time dependence to the 
persistence time distributions, even though the stress re- 
laxation depends on both times and remains finite for 
the duration of the experiment. This is in direct contra- 
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FIG. 6. Relaxation dynamics during the strain step experiment, (a) Strain modulus G = a/e versus time for e = 0.02 (dashed) 
and 0.04 (dotted). Distributions of (b) first hop times, (c) persistence times, and (d) displacements for the undeformed sample 
(solid lines) and e = 0.02 (dashed) and 0.04 (dotted). In all plots t w = 750 (black), 7500 (red), and 75000 (green). Straight 
lines indicate power law with the given slopes. 



diction to the Eyring model of mechanically accelerated 
dynamics, which says that the stress decreases transition 
barriers and hence increases mobility [l|. 



E. Acceleration Factor 

In this section, we develop a method to evaluate the 
complex transformations of the relaxation time distribu- 
tions that were observed for the three protocols, and to 
determine which deformation variables best describe the 
accelerated dynamics. We assume that the deformed re- 
laxation times td are a function of the undeformed relax- 
ation times t u as well as some undetermined parameters 
p which could include, for instance, the stress, strain, 
and/or strain rate: td — td(t u ,p). If the distribution 
of hop times in the undeformed case is p u {t u ), then the 



distribution in the deformed sample becomes: 

Pd(t d (t u ,p)) = Pu(tu)^r (2) 
aid 

Integrating on both sides, we obtain 

[ d p d (t)dt= f Pu (t)dt (3) 
Jo Jo 

where the integrals in Eq. ([3]) are just the cumulative 
distributions P u (t u ) and Pd{td) of the first hop times. 
The acceleration factor can therefore be defined as t u /td, 
where t u and t d are the times when the undeformed and 
deformed cumulative distributions are equal: Pd{td) — 
Pu(t u )- In this way, we obtain the acceleration of the 
entire spectrum of relaxation times during the experi- 
ment and we can determine the form of t d (t ui p), and 
identify the relevant parameters p. This analysis can be 
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FIG. 7. Cumulative probability distribution of the first re- 
laxation event for the undeformed glass (solid line) and un- 
der deformation with a step stress a = 0.3 (o), a step strain 
e = 0.01 (A), and a constant strain rate e = 8.5 x 10 -7 (□) 
for t w — 75000. Dashed line shows the cumulative for the 
undeformed t w = 750 glass. 

performed on both the first hop time and the persistence 
time distributions. 

The cumulatives of the first hop times for represen- 
tatives of the three deformation protocols are shown in 
Fig. [7] It is clear from this figure that the accelera- 
tion factor is unique to each deformation protocol. The 
step strain seems to cause a simple shift in the cumu- 
lative distribution with respect to the undeformed case, 
which means that all relaxation times are accelerated by 
a constant factor. There is also experimental evidence 
for a simple rescaling of the relaxation times under a 
step strain. O'Connell et al. found that the stress relax- 
ation after different step strains could be superimposed 
by shifting the curves in time, that is, the glasses obeyed 
"strain-time superposition" 13]. Alternatively, in both 
the stress step and the constant strain rate experiments, 
the cumulatives for the deformed samples continue to di- 
verge from the reference curve at longer times. This is 
just another way to see that the distribution of relaxation 
times is narrowed: long relaxation times are accelerated 
more than short ones. 

We have already seen that the global stress is not 
a good variable to describe the dynamical acceleration, 
and the fact that the step strain causes a simple shift in 
P(ti) suggests the strain as a promising alternative. In 
Fig. [51 the acceleration factor computed from the first 
hop time distribution i" / tf is plotted parametrically ver- 
sus the strain e for all protocols. The acceleration factor 
for the three different protocols and several wait times 
collapse quite well onto a common curve. No similar 
collapse could be found for the acceleration factor versus 
any other macroscopic deformation variable. The acceler- 
ation factors computed using representative simulations 
with chains 100 beads long are also shown in Fig. [5] and 
have exactly the same form as the short chain results. 
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FIG. 8. Acceleration factor t u /td as a function of global strain 
e for three different deformation protocols. Stress step: a = 
0.3 (o), 0.4 (□), 0.5 (A); constant strain rate: e = 8.9 x 10~ 6 
(0), 8.9 x 10~ 7 (V); strain step (<). For each: t w = 75000 
(black), 22500 (red). Also shown is the acceleration factor for 
the persistence times for constant strain rate deformation at 
e = 8.9 x 10~ 7 (blue *). Solid symbols show results for chains 
of length 100. 



A strain-controlled acceleration factor also explains the 
observed changes to the persistence time distributions. 
We showed in Fig. [T] that individual particle trajectories 
exhibit very little non-affine displacement except during 
a hop. This suggests that the local strain tracks with 
the global strain in between relaxation events, and that 
the relevant parameter in analyzing the effect of deforma- 
tion on the persistence times is not the total strain since 
the onset of deformation, but the accumulated strain in 
between subsequent hops. With this interpretation, we 
understand that the strain step causes no change to p(r) 
because there is no further strain after the initial step. 
Note that this does not imply that there are no residual 
effects of the step strain on the dynamics after the first 
hop has relaxed the local strain: the persistence time dis- 
tribution p(r) is narrower than the first hop time distri- 
bution in an aging glass, therefore particles displaced by 
the deformation are more likely to hop again than those 
that have not. In the constant strain rate experiment, 
we can compute the acceleration factor of the persistence 
times T u /rd from the cumulatives of p(r) in the same way 
as before. Plotted versus the global strain increment in 
between hops (in this case, e is proportional to the persis- 
tence time), Fig. [8] shows that the acceleration function 
for the persistence times is identical to that of the first 
hop times. Although it is more difficult to quantify, the 
dependence of p(r) on the wait time and the measure- 
ment time in the step stress experiment can be explained 
by the fact that the strain explicitly depends on these 
same variables. 

We note here that the acceleration ratios of only the 
longest wait time curves could be evaluated with great 
accuracy. Our method of identifying hops, whereby the 
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average position and standard deviation are calculated 
over a time window, leads to a minimum time resolution 
where individual hops can be identified. The relaxation 
time distributions p{t\) and p(r) are not affected by this 
resolution limit, as discussed in more detail in the Ap- 
pendix. However, the fact that p(ti) and p(r) cannot be 
measured at very short times does presents a problem 
in calculating the cumulatives, especially when there are 
many short time hops, such as for shorter wait times or 
very rapid deformations. In this case, because there are 
many hops at short times that are not detected, the cu- 
mulative increases abruptly from zero at the beginning 
of the measurement, rather than showing the smooth 
increase typical of the cage escape regime (see Fig. [3 
dashed line). The curve is effectively shifted downward 
from where it should be because the integrated effect of 
the shortest time hops is not present. The acceleration 
factor is thus most accurate at long wait times, where 
the cumulative increases slowly from zero, and where 
the short time dynamics in particular are not greatly ac- 
celerated, as in the constant strain rate experiment. In 
the constant strain rate experiment, the acceleration fac- 
tor for all wait times can be evaluated, and the curves 
are found to collapse. In the step strain experiment at 
very large strains, the short time dynamics are acceler- 
ated to such an extent that even the longest wait time 
sample may show truncation effects. To obtain a bet- 
ter estimate of the acceleration factor in this case, we 
additionally used curve fitting to find the factor a where 
Pd(axti)/a most closely resembled p u (ti) and vice versa, 
fitting axp u (ti/a) to Pd{t\). The short time resolution of 
our technique may also be the cause of the widening dis- 
placement distributions with strain rate. If several hops 
occur within our minimum resolution, they are detected 
as a single hop, possibly with a wider displacement. It is 
unclear from our results to what extent the elementary 
hop displacements actually grow with the strain rate. 

The hop statistics provide a powerful way of evaluat- 
ing the full relaxation time distributions and their history 
dependence. Unfortunately, the first hop time and persis- 
tence time distributions are probably not directly acces- 
sible in real polymer glasses. However, it may be possible 
to calculate the acceleration factor from the autocorre- 
lation functions more typically used to measure glassy 
dynamics. In Fig. [§1 we plot the intermediate scattering 
function C(t,t w ) — (exp[— iq • (r(t, t w ) — r(0,t w ))}) un- 
der the same deformation conditions as the cumulatives 
in Fig. [7J Once again the affine portion of the deforma- 
tion is removed and, to minimize the effect of possible 
strain hardening, C is calculated only in directions per- 
pendicular to the direction of strain. It can be seen that 
these curves are qualitatively similar to 1 — P(ti). This 
close connection is to be expected: 1 — P(ii) measures 
the number of particles that have not yet hopped by time 
ti, which is another measure of the particle autocorrela- 
tion. The intermediate scattering function additionally 
includes the effects of multiple hops, since each hop need 
not be very large. This effect causes a steeper decay in 
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FIG. 9. Incoherent scattering function C(t,t w ) with |q| = 
7.0 averaged over the non-deformed (y and z) directions for 
the undeformed glass (solid line) and under deformation with 
a step stress a — 0.3 (o), a step strain e = 0.01 (A), and 
a constant strain rate e = 8.5 x 10~ 7 (□) for t w = 75000. 
Dashed line shows C(t, t w ) for the undeformed t w — 750 glass. 
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FIG. 10. Acceleration factor calculated from the incoherent 
scattering function as a function of global strain e for three 
different deformation protocols. Stress step: a = 0.3 (o), 
0.4 (□), 0.5 (A); constant strain rate: e = 8.9 x 10" 6 (0), 
8.9 x 10" 7 (V); strain step (<). For each: t w = 75000 (black), 
7500 (red). 



C(t,t w ) than in 1 — P(ti), especially when the subse- 
quent hops are also accelerated. An acceleration factor 
t u /t d can be defined for the scattering function decay 
in the same way as for the cumulatives through the re- 
lationship C u (t u ) = Cd(t d ), where C u and Cd are the 
undeformed and deformed correlation functions respec- 
tively. This acceleration factor is shown in Fig. 1 1 01 and 
also approximately collapses as a function of the strain. 
The curves increase nearly exponentially for all strains, 
and more rapidly than the acceleration function calcu- 
lated from the first hop time distributions. The scatter 
in this plot is significantly higher than for the accelera- 
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tion factor calculated from the cumulative distributions, 
for reasons that are not entirely understood. However, 
we believe this measure holds promise for confirming ex- 
perimentally whether the strain is indeed a good way 
to describe deformation-induced acceleration of the re- 
laxation dynamics in glasses. Our analysis evaluates the 
full shape change of the particle auto-correlation function 
due to deformation, rather than reducing the data to an 
average mobility and width of the correlation decay. 



V. DISCUSSION AND CONCLUSIONS 

The statistics of single particle hopping dynamics have 
been evaluated under three different deformation proto- 
cols. All three deformation modes accelerate the segmen- 
tal dynamics, although the exact transformation of the 
relaxation time distributions is specific to the mode of 
deformation. The universal collapse observed in Fig. [5] 
of the acceleration factor as a function of the strain is 
strong evidence that the local strain is a good variable 
to describe the influence of deformation on structural re- 
laxations, rather than the global stress as postulated by 
the Eyring model. A compelling reason for this is sug- 
gested by the behaviour of the stress and strain at a lo- 
cal level. While the local strains track very well with 
the global strain in between relaxation events, the global 
stress does not predict the evolution of stress on a coarse- 
grained scale [3J] . The local distributions of moduli and 
stress surely account for the wide distribution of relax- 
ation times 33] , however the strain appears to be a better 
predictor of how these relaxation times are accelerated. 

By tracking the individual particle relaxations, our 
analysis permits direct observation of the narrowing of 
the relaxation time spectrum during creep and constant 
strain rate deformation, which previously was inferred 
only indirectly from stretched exponential fits to corre- 
lation functions [l8j]. The strain-controlled acceleration 
factor also provides a simple, intuitive explanation for 
the narrowing: particles with long relaxation times ex- 
perience more acceleration than particles that hop sooner 
simply because they experience more strain in the long 
period of time in between relaxations. Alternatively, a 
step strain leads to a simple shift in the first hop times 
(no narrowing), and no change at all to the persistence 
times. The fact that the persistence times are unchanged 
in the step strain protocol suggests that the strain at the 
microscopic level is released after a relaxation event, or 
reset to zero. This interpretation is further confirmed by 
the fact that the acceleration of the persistence times for 
the constant strain rate experiment is identical to that 
found for the first hop times if, instead of the total strain, 
the strain in between subsequent hops is taken to be the 
relevant mechanical parameter. 

Our results do not imply that the average segmental 
dynamics in the glass should be monotonically increased 
with the strain. In their experiments measuring the seg- 
mental dynamics during creep at stresses above yield, 



Lee et al. found that the average segmental mobility in- 
creases with strain at first, becoming a maximum where 
the strain rate was largest, and then decreases again in 
the strain hardening regime (l8j . This is completely con- 
sistent with the strain-controlled acceleration factor pre- 
sented here. The acceleration factor increases with the 
strain only in between subsequent hops; as discussed pre- 
viously, after a relaxation event the local "strain clock" is 
reset to zero. As a result, at the beginning of a deforma- 
tion experiment, the local strains track closely with the 
global strain; however, after some time, the local strains 
in the glass become spatially heterogeneous, being higher 
in regions that have not yet relaxed, and lower in regions 
that have recently experienced a relaxation. In the steady 
state flow regime, our simulations show that every par- 
ticle has hopped at least once and usually many times. 
The advantage of molecular dynamics simulations is that 
the individual trajectories reveal exactly how much time 
has elapsed between local relaxations and thus how much 
strain has accumulated in each local region. 

It is interesting to note that the strain was chosen 
as control variable in the Soft Glassy Rheology (SGR) 
model Q. Here the acceleration factor for barrier cross- 
ings increases proportional to exp[fc/ 2 /2x], where k is the 
elastic constant of a mesoscopic element, I a measure 
of local strain on that element and x denotes a noise 
temperature. The local strains increase in tandem with 
the global strain in between barrier crossings, and after 
a hop, the mesoscopic element returns to a position of 
zero local strain. The main elements of the SGR model 
closely resemble our findings for the single particle re- 
laxation times. The measured form of the acceleration 
factor is, however, quite different. It would be interest- 
ing to see if using the measured form of the acceleration 
factor could account for some of the limitations identified 
in this model [301 ] . 

The hop statistics under deformation also present an 
explanation for mechanical rejuvenation. In contrast to 
the case of sub-yield deformations where the response 
continues to show explicit history dependence, the ef- 
fects of aging are lost once the glass has yielded. Non- 
stationary dynamics occurs in glasses because of the wide 
distribution of persistence times. An infinite (r) means 
that an equilibrium state is impossible, and the dynamics 
experiences aging. During yield in both the step stress 
and the constant strain rate protocols, the persistence 
time distributions are significantly narrowed, and sta- 
tionary dynamics are restored. The resulting rejuvenated 
state is, however, not the same as a just-quenched glass. 
The relaxation time distributions depend explicitly on 
the deformation protocol. Further investigation of the 
relaxation time distributions after deformation may shed 
more light on the nature of the rejuvenated state. 

The hop dynamics under deformation also provide 
some information about strain hardening in polymer 
glasses. We have tested two different chain lengths and 
find that the hop rate versus strain rate and the acceler- 
ation versus strain are identical for both chain lengths. 
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Even though strain hardening is significantly more pro- 
nounced in the long-chain glass, the deformation appears 
to affect the hop dynamics in the same way. Strain hard- 
ening is most evident in the displacement distribution, 
which becomes anisotropic and significantly narrower in 
the strain direction. Further study of the hop dynamics 
during strain hardening may provide further insight into 
this very important behaviour of polymer glasses. 

In conclusion, by detecting individual segmental re- 
laxations, the full transformation of the relaxation time 
distributions under deformation could be quantified and 
correlated with the global deformation parameters in a 
model polymer glass. We find that the accelerated seg- 
mental dynamics are fully explained by the global strain, 



and we present a possible experimental test of these find- 
ings though the definition of the acceleration factor for 
any segmental autocorrelation function such as the inco- 
herent scattering function. 
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Appendix: Identifying hopping events 
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FIG. 11. (a) A typical particle trajectory in the glass. Red 
curve shows raw position data in the x-direction, and the 
black curve shows the running average, (b) Standard devi- 
ation in the three dimensional particle position over the av- 
eraging window. Hops are identified by a threshold in the 
standard deviation, shown here as a horizontal dashed line, 
and marked in both frames as vertical dashed lines. In panel 
(a) the first hop time, persistence times, and particle displace- 
ments are also labeled. 



Particle positions are written to file with a period of 
tdump- The trajectories show rapid vibrations about a 
mean position, punctuated by the occasional hop to a 
new location. To detect the hopping transitions, the 
average and standard deviation of each particle posi- 
tion is computed over N avg snapshots at times ti — 
N aV g x tdump X (i — 1/2). The trajectories and stan- 
dard deviation are shown in Fig. [IT] The value of N avg 
used for determining the mean and standard deviation is 
on the order of 40 — 60 snapshots for the work presented 
here. Increasing N avg by a factor of 10 (from 40 to 400), 
causes negligible changes to the hop statistics but seri- 
ously decreases the time resolution of the distributions. 
The criteria for finding a hop is that the standard de- 
viation of the particle position over N avg snapshots is 
greater than a threshold: a > Uth- To choose a suitable 
threshold, we look at the distribution of the standard de- 
viations over the full time series of all of the particles in 
the system (see Fig. [T2|) . There are clearly two regimes 




FIG. 12. Probability distribution of the standard deviation a 
in the particle position over an averaging time window of 400 
for all particle trajectories. Dashed line shows the threshold 
ath for detecting a hop. 



of motion in p(a): a narrow vibrational regime and a 
hopping regime. tj t h is chosen to be at the shoulder in 
between these regimes and is shown in Fig.[12]as a dashed 
line. The effects of the threshold on the distributions is 
discussed below. 

The times at which a hop occurs are determined by the 
centre of the standard deviation peak, and the displace- 
ments are found from the positions on either side of the 
peak. The distribution of displacements is easily com- 
puted by binning the data to form a histogram. Because 
the relaxation time distributions are very broad, the bin 
sizes for p(t\) and p(r) increase logarithmically with time 
in order to obtain good statistics. The finite timespan of 
the simulation means that this distribution is necessarily 
truncated. In other words, some of the particles will not 
have hopped at all during the experiment. For the first 
hop times t\ , this can be accounted for by normalizing to 
the total number of particles in the ensemble, rather than 
the number of first hops detected. The persistence time 
r is defined as the time in between all subsequent hops. 
Accounting for the effects of truncation on this distribu- 
tion is somewhat more complex than in the case of p{t\). 
If all of the hops are assembled into a single distribution 
function as described above, there is a sharp downturn 
in p(r) near the maximum time observed t max . The dis- 
tribution of persistence times that we observe p b s (t) is 
weighted by the probability that we will observe it within 
t < t max , which becomes low for cages originating near 
the end of the simulation time. This effect can be ac- 
counted for using multiple convolutions of the relaxation 
time distributions; however since the downturn is very 
sharp, it is in practice easier to simply discard the very 
long time data. 

The effect of the dump frequency on the distributions 
was investigated by recording the same data at three dif- 
ferent values of tdump- Figure [JJJ] shows the effect of tdump 
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FIG. 13. The distributions (a) p(ti), (b) p(r) and (c) p(dx) for tdump = 3.75 (circle), 37.5 (square) and 375 (triangle). 




FIG. 14. The distributions (a) p(ti), (b) p(r) and (c) p(dx) for oth = 0.15 (square), 0.25 (circle), and 0.35 (triangle). Solid 
lines in (c) indicate power laws with slope of -1.5 and -1.1. 




FIG. 15. (a) The distribution of hop correlations, (b) p(r), and (c) p(dx) computed with all of the hops detected (circles), and 
with only uncorrelated hops (squares). 



15 



on (a) p{h), (b) p(r) and (c) p(dx). Changing t dump does 
not affect the first hop time distribution p{ti), aside from 
the obvious relationship between tdump and the minimum 
t\ that can be observed. Distributions recorded at differ- 
ent values of tdump are in perfect agreement where their 
time spans overlap. The persistence time distribution 
p(r) has the same slope for all values of tdump-, how- 
ever the curves differ by a normalization factor which de- 
pends on the minimum time resolution in r. An absolute 
normalization factor may be obtained by using a similar 
method to that used to normalize the distributions p(ti): 
the time between the first and second hops are found for 
a subset of N particles where t\ is very short. Not all of 
these particles will experience a second hop, and if the 
distribution is normalized to N rather than the number 
of second hops, a normalized distribution is obtained that 
does not depend on tdump- This method requires throw- 
ing away most of the persistence time data, however, and 
is therefore only used to properly normalize the full dis- 
tribution which includes all measured persistence times. 
The displacement distribution p{dx) becomes somewhat 
wider with increased tdump due to the effects of unde- 
tected hops on timescales shorter than tdump- The dis- 
tributions compared here have identical tdump- 

The effect of the threshold on the standard devia- 
tion was also investigated using three different values: 
Oth = 0.15, 0.25, and 0.35. The distributions computed 
using these thresholds are compared in Figure [T?] Chang- 
ing the threshold results in quantitative, rather than 
qualitative changes to the distributions. In all cases, the 
hop time distributions take the form of power laws, al- 
though the exponent is modified by the threshold (from 



about -1.1 to -1.5 for p(r)), and the first hop time distri- 
bution depends on t w whereas the persistence time distri- 
bution does not. However, if the threshold is very small, 
the influence of wait time is decreased because of noise 
(vibrations rather than hops) that does not depend on 
t w , if 

<j t h is very large, the statistics become poor and 
many legitimate hops are discarded. The hop displace- 
ments predictably become larger with increasing thresh- 
old. Using a continuous time random walk |46| . we find 
that the mean squared displacement and van Hove func- 
tion are self-consistcntly described by the hop statistics 
obtained using a range of thresholds near the shoulder of 
p(cr) (Fig.[T2l. 

In order to model the trajectories as a continuous time 
random walk we assume that the hops are completely 
uncorrelated. To ensure that this is a reasonable approx- 
imation, the correlations between adjacent hops in the 
trajectories are evaluated. A correlation function can be 
defined for each pair of displacement vectors 

Q = llr^lllNr (A ' 1} 

The distribution of these correlations are evaluated for 
all pairs of hops, as shown in Fig. [TSJa). A completely 
uncorrelated random walk would have a flat distribution; 
however, we see that there are a number of backward cor- 
related (C = —1) and forward correlated (C — 1) hops. 
If these hops are removed from the time series data and 
the distribution of persistence times and displacements 
are computed with only the uncorrelated hops, we see 
no appreciable difference in the hop statistics (Fig. ITBTb') 
and (c)). 



